#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<gmp.h>
#include<pthread.h>
#include<time.h>

//algoritmo quartic convergence(1984)
//Then pk converges quartically to π; that is,
//each iteration approximately quadruples the number of correct digits.
//The algorithm is not self-correcting;
//each iteration must be performed with the desired number of correct digits of π.

#define NUMTHR 8

pthread_t t[NUMTHR];/*Cria as threads*/
mpf_t arg[10];/*Cria as variaveis da gmp*/

void * borwein(void *);

int main(void){
	int i = 500, j=0, k=0;
	time_t t1, t2;
	FILE *fp;
		
	fp = fopen("pi_input.txt", "r");
	if(fp == NULL){
		printf("Erro na abertura do arquivo!\n");
		return -1;
	}

 	time(&t1); /*pega o tempo nesta etapa da execucao para posterior comparacao*/
    /*Inicializa as variaveis da gmp*/
       for(j=0; j<10;j++){
       	mpf_init(arg[j]);
       }

    /*define a precisao*/
       mpf_set_default_prec(33219280);

	/*atribui valor 0 as variaveis gmp*/
       mpf_set_ui(arg[0],0);
       mpf_set_ui(arg[1],0);
       mpf_set_ui(arg[2],0);
       mpf_set_ui(arg[3],0);
       mpf_set_ui(arg[4],0);
       mpf_set_ui(arg[5],0);

    /*atribui valor da constante pi na variavel gmp*/
	printf("Retorno do input file = %d\n", mpf_inp_str (arg[9], fp, 10)); 
	printf("\nPi = ");mpf_out_str (stdout, 10, 10, arg[9]); printf ("\n\n");
	
    //arg[0] = sqrt(2);
    mpf_sqrt_ui (arg[0], 2);
    //arg[1] = 0;
    //arg[2] = 2 + sqrt(2);
    mpf_sqrt_ui (arg[2], 2);
    mpf_add_ui (arg[2], arg[2], 2);

	while(mpf_cmp(arg[9],arg[5]) <= 0 || i > 0)
	{
  	  for(k = 0; k< NUMTHR; k++ ){
		    pthread_create(&(t[i]), NULL, borwein, arg);
		    pthread_join(t[i], NULL);
	  }
	  i--;
	}    
	time(&t2);
		  printf("Iteracoes = %d, tempo calculado = %ld\n", (500-i), (int) t2-t1);
		  gmp_printf("Pi %.20Ff \n",arg[5]);    
		

    return 0;
}

void * borwein(void *argument){
	mpf_t *arg = (mpf_t*) argument;
	
        //an = (sqrt(arg[0]) + (1/sqrt(arg[0])))/2;
        mpf_sqrt(arg[6],arg[0]);
        mpf_ui_div (arg[7], 1, arg[6]);
        mpf_add (arg[8], arg[6], arg[7]);
        mpf_div_ui (arg[3], arg[8], 2);

        //bn = ((1+arg[1])*sqrt(arg[0]))/(arg[0]+arg[1]);
        mpf_add_ui (arg[7],arg[1], 1);
        mpf_mul (arg[8], arg[7], arg[6]);
        mpf_add (arg[6], arg[0], arg[1]);
        mpf_div (arg[4], arg[8], arg[6]);

        //pn = ((1+an)*p0*bn)/(1+bn);
        mpf_add_ui (arg[6], arg[3], 1);
        mpf_mul (arg[7], arg[6], arg[2]);
        mpf_mul (arg[8], arg[7], arg[4]);
        mpf_add_ui (arg[6], arg[4], 1);
        mpf_div (arg[5], arg[8], arg[6]);

		//arg[0] = an;
		mpf_set (arg[0], arg[3]);
		//arg[1] = bn;
		mpf_set (arg[1], arg[4]);
		//p0 = pn;
		mpf_set (arg[2], arg[5]);
    
}
